A nonlinear equation for ionic diffusion in a strong binary electrolyte 
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The problem of the one dimensional electro-diffusion of ions in a strong binary electrolyte is 
considered. In such a system the solute dissociates completely into two species of ions with unlike 
charges. The mathematical description consists of a diffusion equation for each species augmented 
by transport due to a self consistent electrostatic field determined by the Poisson equation. This 
mathematical framework also describes other important problems in physics such as electron and 
hole diffusion across semi-conductor junctions and the diffusion of ions in plasmas. If concentrations 
do not vary appreciably over distances of the order of the Debye length, the Poisson equation can 
be replaced by the condition of local charge neutrality first introduced by Planck. It can then be 
shown that both species diffuse at the same rate with a common diffusivity that is intermediate 
between that of the slow and fast species (ambipolar diffusion). Here we derive a more general 
theory by exploiting the ratio of Debye length to a characteristic length scale as a small asymptotic 
parameter. It is shown that the concentration of either species may be described by a nonlinear 
integro-differential equation which replaces the classical linear equation for ambipolar diffusion but 
reduces to it in the appropriate limit. Through numerical integration of the full set of equations it 
is shown that this nonlinear equation provides a better approximation to the exact solution than 
the linear equation it replaces. 
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PACS numbers: 82.45.Gj,82.45.-h,82.70.y,47.57.eb,47.57.jd,72.20.Dp 



The one dimensional electro-diffusion equations de- 
scribing the evolution of ionic concentrations n(a:,i), 
p(a;, t) in a fully dissociated binary electrolyte with a self 
consistent electrostatic potential <^(a;, t) may be put in 
the form [l|] 



dtn 


- dxindxcj)) = 


dxxn 


(1) 


-'dtp - 


zdx{pdj:(l)) = 


dxxP 


(2) 




dxx(t) = 


—n — zp 


(3) 



where u = Up/un is the ratio of mobilities of species p 
and n and z = Zp/zn is the ratio of their valences. We will 
take n as the species with lesser mobility, so that u > 1. 
Since the two species are oppositely charged, z is neg- 
ative. The above equations have been rendered dimen- 
sionless by scaling the concentrations by the concentra- 
tion of n at infinitely far points (n*) so that n{oo, t) = 1, 
the co-ordinate a; by a length scale A* related to De- 
bye length, the time by a diffusion time A^/-D„ and the 
potential cf) by the thermal energy (kBT)/(eZn)- In the 
above, Z3„ and Dp are the diffusivities of the two species 
which are related to the mobilities Un,Up through the 
Einstein relation Z)„/u„ = Dp/up = ksT, where /cb is 
the Boltzmann constant and T the absolute tempera- 
ture. The length scale A* referred to above is defined by 
the expression A^ = (eonkBT) / (z^e^rii,) in terms of the 
electronic charge (e), permittivity of vacuum (eq) and di- 
electric constant (k). It is related to the Debye length in 
a homogeneous charge neutral electrolyte where the less 
mobile species has a dim ensionless concentration of n as 
AD(n) = A*/^n(l - z). 

Eq.([T])-(l3]) and its modifications also describe other 
transport problems of interest in various areas of physics. 
For example, if 1/1(2;, t) = —EqX + (J)' where Eq is an exter- 



nal field and (j)' is the perturbation in the potential and if 
a term —N[x) representing the density of free charges is 
added to the right hand side of Eq. ([3]) , then these equa- 
tions become the Van Roosbroeck equations describing 
the migration of charge carriers in solid state physics, 
where n, p and N represent the concentrations of elec- 
trons, holes and dopants [2]. If N{x) is interpreted as the 
density of ion exchange sites within a membrane, then 
the same set of equations also describes various filtration 
processes such as electrodialysis [1] . In plasma physics Q 
one often encounters situations where the electrons and 
positive ions may be regarded as two interpenetrating 
charged fluids of different diffusivities coupled by a self- 
consistent electric field. It is then described mathemat- 
ically by Eq. (HI) -([3]). The generalization of these equa- 
tions to three or more species in the presence of an ap- 
plied electric field Eq describe the separation, of charged 
macromolecules in capillary electrophoresis |^ and other 
separation processes based on electric charge. 

One of the simplest experimental realizations of 
Eq. (HI)-® is the "Liquid Junction" problem, where a 
barrier initially separates two semi-infinite regions of a 
binary electrolyte (such as sodium chloride in water) of 
different concentrations, and subsequently, at time i = 0, 
the barrier is removed. Since generally the positive and 
negative components would diffuse at different rates a 
polarization vector and consequently a measurable "Liq- 
uid Junction Potential" (LJP) appears across the inter- 
face. Planck Q derived an expression for the LJP on 
the basis of the local charge neutrality assumption, an 
idea that has since become a central paradigm in many 
problems involving electro-diffusion. The physical basis 
of the approximation is that the electrostatic force that 



would arise if a charge separation was realized is so strong 
that in practice it precludes any departure from electro- 
neutrality anywhere in space. 

Mathematically, the local electro-neutrality assump- 
tion amounts to neglecting the term Oxxfj) on the left 
hand side of Eq. ([3]), so that it reduces to p — ~n/z. 
This relation can now be used to eliminate the terms in- 
volving (j) in Eq. ([T]) and ([2]). It then follows that n (as 
well as p = —n/z) obeys the diffusion equation 



dtn^ Ddxx 



(4) 



where D = u{l — z)/(l — uz), a result due to Hender- 
son fo, '?']. Since z <0, clearly u> D > 1. The coulomb 
attraction between the two species enhances the rate of 
spreading of the slower species and reduces that of the 
faster species so that both diffuse at an equal, ambipolar 
(that is, the same for either charge) rate. In the anal- 
ysis outlined above, though the term in (p is dropped 
in Eq. ([3|), it must be retained in Eq. ([T])-(l2]), a situ- 
ation that appeared self contradictory and led to some 
controversy until Hickman [8|, |9| provided a proper inter- 
pretation within the framework of an asymptotic theory 
based on expansion in the small parameter /cAu, where 
fc is a characteristic wave number and Xd is the Debye 
length. Since the Debye length Xo is normally a very 
small quantity (about 3 nm in a 0.01 M solution of a 
common salt like sodium or potassium chloride) in most 
applications the charge neutrality assumption is very ac- 
curate. However, it could be violated in many modern 
applications such as in nanochannels where one or more 
geometric dimensions may be of the order of nanometers. 
Such departures from charge neutrality may give rise to 
new effects. It would therefore seem worthwhile to first 
study these effects in the context of the simple model 
system represented by Eq. ([J)-©. 

In this paper we use the method of multiple scales [l3| 
to reduce Eq. (IT])-(|31) to a nonlinear one dimensional sys- 
tem in the limit of long (compared to the Debye length) 
length scales and slow (compared to the diffusion time 
over a Debye length) time scales. We show that at the 
lowest order the equation for ambipolar diffusion is re- 
covered. If the asymptotic theory is continued to the 
next order, a nonlinear integro-differential equation for 
the concentration n emerges. Numerical solutions of this 
reduced equation is seen to agree better with that of the 
full system of Eq. ([I])-© and lead to effects not captured 
in the lowest order theory based on local charge neutral- 
ity. 

We would like to study the behavior of Eq. (HJ-Q at 
long length and time scales, under the boundary condi- 
tions that n,p, (j) all approach constant values and respect 
local charge neutrality as a; —J' ±oo. Thus, we are consid- 
ering passive diffusion in the absence of imposed electric 
fields and currents. Following the usual procedure [10|, 
we introduce slow variables ^ — \fix and r = ei and 
suppose that n,p^(f> depend solely on ^ and r. Then 



Eq. dJ)-© reduce to 

drU — d^{nd^(j)) — d^^n (5) 

u~^drp— zd^{pd^(j)) = d^^p (6) 

e9jj0 = —n — zp (7) 

where e is a small parameter. Expanding all dependent 
variables in asymptotic series in e, such as n = nQ + eni + 
e^n2 + • • • , substituting in Eq. ©-((T]) and equating the 
coefficients of e on both sides, we get a series of equations 
starting with the lowest order ones: 

driiQ - d^{nQd^(f)f)) = %no (8) 

u~^drPo - zd^{pod^<j>Q) = d^^po (9) 

-no - zpo ^ (10) 

If the last equation is used to eliminate the second terms 
from the first two, we get 



9rno = Dd^^no 

Po = -no/z 



(11) 
(12) 



with D = u(l — z)/[l — uz), representing ambipolar dif- 
fusion and if the time derivative terms are eliminated 
instead and the result integrated, we get 



u- 1 
1 — uz 



In no. 



(13) 



If the above equation is evaluated at a: = — oo, we get 
Planck's formula 5] for the potential drop across a liquid 
junction. At the next order in e, we have 

9r"-i -9^("i950o)-9^(noa40i) = %ni (14) 
u^^drPi- zd^{pid^(j>o)~ zd^{pod^<j)i) = %pi (15) 

(16) 



d> 



'^^(Po 



-ni - zpi 



If we add Eq. ^^ and ([T5]) and use Eq.([T0l), the term 
in 01 is eliminated. We then obtain after substituting c/jq 
from Eq. (fT5|) and after some algebra; 



drHi = Dd^^ni + ad^^^^{\iino) - (3d^^{d^\nno)'^ (17) 

u — 1 

%(ln?io), (18) 



ni 
Pi = 



z z(l — uz) 
where a and /3 are positive constants defined by 
uz{u — 1)^ 



/3 



(1-^2)3 

u{u — 1)(2 — z ~ uz) 
2(1-Mz)3 



(19) 
(20) 



If we add Eq. (|TT|) to e times Eq. (ITTI) and transform back 
to our original variables x and t, we get 

dtn = Ddxxn + adxxxxi'^'^i^n) - I3dxx{dx Inn)^ (21) 

with an error which is higher than order e^. In arriv- 
ing at Eq. (PT|) we have replaced the terms ad^^^^{\nno) 



by ad^^^^{lnn) and — /39^^(9^ Inno)^ by — /3i9^^(9^ Inn)^, 
which is justified because doing so introduces an error in 
Eq. ()2ip that is higher than order e^. Similarly, combin- 
ing Eq. ^ with Eq. ^ we get 



1 



P 



z z(l — uz) 



dxxi^^n) 



(22) 



with an error of order e^. The last term in Eq. ip^ is due 
to the departure from charge neutrality. 

Let us now consider some of the consequences of 
Eq. (|2T|) and ([22l) . First, we consider weak perturba- 
tions from the constant values at infinity: n = 1 + n' 
and p = —1/z + p' , where \n'\, \p'\ <C 1. This gives a 
hyper-diffusion equation for n': 



dtti' = Ddxxn' + adxxxxn' ■ 



(23) 



Therefore, solutions of the form exp[ai-f ifca;] have growth 

Since a > 0, high wavenumber 



rates a — —k^D 
modes k > fc^ax 



k" 



y/{D/a) are unstable. We will show 
that this instability is entirely spurious. Indeed, Eq. (PT|) 



is not even valid for modes k > k„ 



1 since these solu- 



tions violate the premise fc <C 1. The instability arises as 
a consequence of truncating the asymptotic series, as can 
be seen by considering the linearized version of Eq. p^- 



© with u-^ 


= 0: 










dtn 


- dxx4> = 


dxx-n', 


(24) 






dxx4> = 


dxxP', 


(25) 






dxx4> = 


—n' — zp' . 


(26) 



If Eq. (|25)) is integrated and the result substituted in 
Eq. (|26p an exact expression for may be found in terms 
ofn': 



-{l + z'^OxxY^n' 



+ 00 



G{\x-y\■,^/^)n'{y)dy, 



(27) 
(28) 



[1 - z-^dxx + z-'^dxxxx + •••]«' (29) 



where G{x]m) = — cxp(— TO|x|)/(2m) is the Green's 
function of the Helmholtz operator dxx — 'm? ■ If the ex- 
act inversion, that is Eq. ([28| . is substituted in Eq. (|24|) 
one gets an integro-differential equation which in Fourier- 
space yields a growth rate 



= -e 



1 



-fc'^ ( 1 - - 1 + ^ 



(30) 



The exact formula for a indicated by the — sign above 
shows that there is no high wavenumber instability if the 
exact inversion, Eq. (P5|. is employed. The second part, 
indicated by the ~ sign, shows the result of the approx- 
imate inversion, Ea. (j29p . based on the low wavenumber 
approximation. In this case there is a high wavenumber 
instability if the asymptotic series is truncated after an 
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FIG. 1: The quantity A = (2Dy^da'^/dt - 1 (lower quad- 
rants) where a^ is the variance of the distribution n{x,t), and 
the excess kurtosis 7 (upper quadrants) determined using the 
models (II) [dashed] and (III) [solid] . For model (I) , 7 = A = 
at all t. The initial state is n{x,0) = 1 + Aexp{—x^ /2) and 
parameters are A = 1.0, u = 6, z = —1. 



even number of terms. That is, the root cause of the 
instability is the oscillatory approach to the limit of the 
series indicated on the right hand side of Eq. ([5D|) . 

The above analysis suggests a simple way of modifying 
Eq. (|21[) without violating its asymptotic validity, but at 
the same time eliminating the spurious high wavenumber 
instability. We recognize that the term adxxxx(}'nn) — 
dxx [oidxx In n) in Eq. (j2ip most likely resulted from a 
truncated expansion of the Helmholtz operator: dxx[^ — 
adxx]~^ ~ dxx{l + adxx + ■■■) = dxx + adxxxx + ■ ■■ 
and simply replace adxx in Eq. (j2ip by (1 — adxx)~^ ~ 1- 
Thus, we get the nonlinear integro-differential equation: 



dtn = dxxJ^[n\ 



(31) 



where J^ is the functional: 

^■[ti] = Dn- I3[dxilnn)] 
1 r+°° 



+ 



2^/5 



e-^''-y^/'^\n 



n{y) 
n{x) 



dy. (32) 



Eq. ([5T|) is asymptotically equivalent to Eq. (1^ since 
they differ only by terms of order higher than e^ but it is 
free from the spurious high wavenumber instability. By 
virtue of its construction it also has the property that 



when 



— 0, the linearized version of Eq. ([3T|) is the 



exact solution of Ea. (l^^ - (P5)) . Thus, when u is large 
and amplitudes are small, the solutions of Eq. (|3T|) match 
closely the true solution, irrespective of the validity of 
the assumption of long length scales and slow time scales 
exploited in the asymptotic theory. The formal condi- 
tion for the validity of Eq. ([3T|) is XDn~^dxn <^ I, where 
the bars on top of the variables indicate that they are 
in dimensional form. In dimensionless notation, the con- 
dition of validity reduces to \dxn\/[rv^/^\/l — z] <C 1. It 



should be noted that, Eq. ((3T|) as well as the lower or- 
der effective diffusivity approximation fails in the limit 
of very low concentrations (n — ^ 0) . This is because the 
Debye length becomes infinitely large in this limit inval- 
idating our premise of long length scales compared to 
the local value of the Debye length. Eq. ((3T|) together 
with Eq. p2|) to calculate the concentration of the faster 
diffusing species are our principal results. 

A numerical method employing a fourth order finite 
difference algorithm for spatial derivatives coupled with 
a fourth order Runge-Kutta time stepping scheme with 
fixed grid and step size was used to solve a time de- 
pendent electrodiffusion problem. An initial condition 
n{x,0) = 1 + Aexp{—x'^/2) was used with A = 1.0 and 
the other parameters were set as u = 6 and z = — 1. The 
problem was solved at three levels of approximation: 

(I) Using the ambipolar diffusion equation, Eq.(|4]). 

(II) Using the nonlinear equation, Eq. (|3T|) . 

(III) Using the full electro-diffusion model Eq. ([I])-©- 

It is clear that with (I) the variance ct^ would increase 
linearly with time, but not so in the case of (II). Further, 
a distribution that is initially Gaussian remains so at fu- 
ture times under (I) but not under (II). Therefore, the 
excess kurtosis 7 is expected to remain zero at all times 
in the case of (I) but not in the case of (II). In Figure 1 
the quantities A = {2D)^^d(j^ /dt — 1 and 7 are plotted 
as a function of time (t) for (II) and (III). In case of (I), 
A = 7 = at all i and this case is not shown for clarity. 
It is seen that at sufficiently long times (I), (II) and (III) 
all approach a common asymptotic limit. However, at 
shorter times, (II) is in better accord with (III) than (I) 
is. 

The qualitative nature of the time variation of 7 and 
A may also be understood on the basis of Eq. ([31]) . To 
see this, we use Ea.p9| to expand the integral operator 
and put Eq. (f3T|) in the form 



dtn = Ddxxn+adxxxx0^^n)-(3dxx{dx'^'D-n)'^ + - 



(33) 



where • • • indicate terms involving sixth, eighth, tenth, 
.... order derivatives. If we further restrict ourselves to 
small amplitudes, n = 1 + n' , with \n'\ <C 1, then Inn — 
n' — ■^{n')'^ + • • • , so that the Eq. ([33l) becomes 

dtn' = Ddxxu' +adxxxxn' --dxxxx{n'Y-(3dxx{dxn'f+- ■ 

(34) 



Here the • • • now include the nonlinear terms of higher 
order. If we multiply Eq. (IMl) by x^ and integrate, we 
can generate ordinary differential equations for the mo- 
ments ruk = /_ x''n'{x) dx, starting with dmo/dt — 
dmi/dt = 0. The neglected sixth and higher derivative 
terms do not contribute for fc < 4. Without loss of gen- 
erality, we can assume mi = so that here cr^ = 7712/7710 
and 7 = 17141710/1712 — 3. By combining the equations for 
the second and fourth moments, we derive 



^-=2D-^ {dxn'fdx (35) 

at mo J_oo 

and 



dj 
'dt 



47D 



12a 



(36) 



where • • • indicate the nonlinear terms that we suppress. 
Now, if a = 0, then either 7 = (if it is zero initially) 
or else 7 monotonically decays to zero. When deviations 
from the ambipolar diffusion limit is considered (a > 0), 
we find that 7 initially increases (when a is small enough 
for the second term to dominate), but as a becomes 
larger, the first term eventually dominates resulting in a 
decrease in 7 towards zero. Further, Eq. ([55)) shows, that 
since /3 > the rate of increase of the variance is slightly 
less than 2D but the deficit in the growth rate becomes 
progressively smaller at large times. This is indeed what 
is observed in Figure 1. 

In summary, a prototypical problem in electro- 
diffusion was considered in the limit where the Debye- 
length was non-zero but nevertheless small compared to 
a characteristic scale of spatial variation. It was shown 
that in this limit the concentration can be described by 
a one dimensional nonlinear integro-differential equation 
which reduces to the linear diffusion equation describ- 
ing ambipolar diffusion if all but the leading order terms 
in the ratio of Debye-length to a characteristic spatial 
scale are neglected. Numerical solution shows that the 
nonlinear integro-differential equation provides a better 
approximation to the true solution. The approach de- 
scribed here could be useful for other problems in electro- 
diffusion where the Debye length is small but may not be 
considered negligible in relation to other length scales. 
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